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Abstract 

One  major  issue  in  the  accurate  solution  of  convection-dominated  problems 
by  means  of  high-order  methods  is  the  ability  of  the  solver  to  maintain 
monotonicity.  This  problem  is  critical  for  spectral  elements,  where  Gibbs 
oscillations  may  pollute  the  solution.  However,  typical  filter-based  stabi¬ 
lization  techniques  used  with  spectral  elements  are  not  monotone.  In  this 
paper,  residual-based  stabilization  methods  originally  derived  for  finite  ele¬ 
ments  are  constructed  and  applied  to  high-order  spectral  elements.  In  par¬ 
ticular,  we  show  that  the  use  of  the  Variational  Multiscale  (VMS)  method 
greatly  improves  the  solution  of  the  transport-diffusion  equation  by  reduc¬ 
ing  over-  and  under-shoots,  and  can  be  therefore  considered  an  alternative 
to  filter- based  schemes.  We  also  combine  these  methods  with  discontinuity 
capturing  schemes  (DC)  to  suppress  oscillations  that  may  occur  in  proximity 
of  boundaries  or  internal  layers.  Additional  improvement  in  the  solution  is 
also  obtained  when  a  method  that  we  call  FOS  (for  First-Order  Subcells)  is 
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1.  Introduction 


A  large  number  of  physical  applications  relies  on  the  accurate  solution 
of  the  transport-diffusion  equation 

|  +  £(«)  =  /,  (i) 

where  q  is  the  concentration  of  the  tracer,  C(q)  =  u  •  Vg  —  V  •  {v's/q)  ,  v  >  o 
is  a  diffusion  coefficient,  u  is  a  known  velocity  field,  and  /  is  a  constant 
source  term.  The  solution  of  (1)  should  respect  two  significant  properties: 
(i)  positivity  must  be  preserved,  and  (ii)  smearing  at  internal  and  bound¬ 
ary  layers  should  not  be  excessive.  These  properties  are  extremely  impor¬ 
tant  in  the  context  of  transport  in  the  atmosphere.  Both  limited-area  and 
global  atmospheric  models  for  weather  prediction  need  monotonic  advection 
of  tracers  and  moisture  variables,  otherwise  the  wrong  amount  of  precipita¬ 
tion  would  be  forecasted.  Simple  microphysics  schemes,  such  as  the  Kessler 
parametrization  [1],  require  three  variables  (water  vapor,  cloud  water,  and 
rain) ,  whereas  more  sophisticated  parameterizations  include  additional  vari¬ 
ables  such  as  ice  and  snow  [2] .  Similarly,  climate  models  require  transport  of 
hundreds  of  tracers,  each  representing  a  different  chemical  species.  Regard¬ 
less  of  the  physical  scales  of  the  model,  tracers  must  remain  positive  since 
the  physical  parameterizations  that  govern  sub-grid  scale  processes  such 
as  auto-conversion  and  sedimentation,  implicitly  assume  such  a  condition. 
These  issues  have  been  addressed  for  both  transient  and  stationary  prob¬ 
lems  (See,  e.g.,  [3])  and,  in  the  context  of  finite  element  methods,  so-called 
stabilized  methods  have  been  an  active  topic  of  research  since  their  intro¬ 
duction  in  the  early  1980s  with  the  Streamline-Upwind  method  of  Hughes 
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and  Brooks  [4].  In  this  paper  we  address  the  problem  of  solving  (1)  by 
high-order  spectral  element  methods  (SEM)  without  losing  the  ability  to 
approach  a  monotone  solution  of  the  problem.  Higher-order  accuracy,  in 
fact,  comes  at  the  price  of  aliasing  phenomena  in  the  solution  [5],  but  the 
anti-aliasing  filters  typically  used  to  give  a  stable  spectral  element  solution 
do  not  respect  conditions  (i)  and  (ii)  described  above.  Therefore,  to  achieve 
monotonic  results  with  high-order  spectral  elements,  we  consider  stabiliza¬ 
tion  schemes  originally  devised  for  finite  elements,  and  focus  on  techniques 
that  can  be  derived  directly  from  subgrid  scale  considerations  as  originally 
defined  in  [6]  and  [7]  in  the  context  of  variational  multiscale  methods.  These 
schemes  assure  stability  by  designing  a  diffusion-type  term  that  is  added  to 
the  Galerkin  formulation  of  the  original  problem. 

The  first  stabilized  schemes  based  on  the  addition  of  a  diffusive  stabi¬ 
lization  term  to  the  Galerkin  equation  are  the  Artificial  Viscosity  methods 
(AV)  [8]  and  the  Streamline-  Upwind  method  (SU)  [4],  AV,  as  Hyper-viscosity 
(HV),  is  often  used  in  atmospheric  and  ocean  modeling  due  to  the  property 
of  preserving  the  correct  energy  cascade  in  simulations  that  involve  turbu¬ 
lence.  The  SU  scheme  uses  the  information  in  the  direction  of  the  flow  to 
add  viscosity  only  in  the  streamline  direction.  Both  methods  use  a  constant 
diffusion  coefficient  that  does  not  typically  change  from  element  to  element. 
A  major  improvement  came  by  introducing  the  residual  of  the  governing 
equation  in  the  definition  of  the  stabilization  term.  When  the  computed 
solution  approaches  the  exact  solution,  the  stabilization  term  should  van¬ 
ish.  This  strategy  is  known  as  residual  weighting  and  generates  a  family  of 
stabilization  methods  used  mostly  in  FEM-based  Computational  Fluid  Dy¬ 
namics  (CFD).  These  schemes,  which  are  consistent  in  that  the  stabilization 
terms  goes  to  zero  as  the  numerical  solution  approaches  the  exact  solution, 


4 


are  considered  in  this  paper.  The  most  commonly  used  are  the  Upwind- 
Streamline /Petrov- Galerkin  (SUPG)  and  the  Galerkin/ Least- Squares  (GLS), 
devised  in  1982  [9]  and  1989  [10],  respectively,  as  a  consistent  counterpart 
to  SU.  GLS  was  designed  as  a  generalization  of  SUPG,  but  in  the  limit 
of  pure  advection,  or  for  piece-wise  linear  elements,  the  GLS  and  SUPG 
methods  are  equivalent.  Stability  analysis  for  these  two  methods  is  detailed 
in  [11,  12,  10].  The  Gradient  Galerkin /Least- Squares  [13]  for  advection- 
diffusion  with  a  reaction  term,  or  the  Unusual  Stabilized  Finite  Element 
Method  (USFEM)  [14,  15]  are  a  few  examples.  In  the  framework  of  high  or¬ 
der  methods,  Petrov-Galerkin  stabilization  was  applied  by  Pasquarelli  and 
Quarteroni  [16]  to  stabilize  the  convection-diffusion  equation  with  the  spec¬ 
tral  method.  Canuto  used  bubble  functions  to  address  the  same  issue  [17] 
(See  also  [18,  19]). 

The  analyses  of  Hughes  [6] ,  Hughes  and  Stewart  [20] ,  and  Hughes  et  al. 
[7]  form  the  unifying  theory  of  all  stabilized  finite  element  methods.  Accord¬ 
ing  to  this  theory,  stabilized  methods  are  subgrid  scale  models  where  the  un¬ 
resolved  scales  are  intimately  related  to  the  instabilities  at  the  level  of  the  re¬ 
solved  scales,  and  thus  should  be  used  in  the  construction  of  the  stabilization 
term.  These  schemes  are  known  as  Variational  Multiscale  (VMS)  methods. 
Details  are  given  in  subsection  2.2.2.  VMS  methods  are  all  residual-based 
methods  that  improve  the  stability  properties  of  the  solution,  and  preserve 
the  accuracy  of  the  underlying  numerical  scheme  [21].  However,  Godunov’s 
theorem  [22]  implies  that  the  latter  property  may  be  violated  in  the  prox¬ 
imity  of  discontinuities  or  strong  gradients.  To  the  authors’  knowledge,  the 
only  application  of  VMS  to  spectral  elements  is  the  work  of  Wasberg  et  al. 
[23]  in  the  context  of  large  eddy  simulation. 

Neither  SUPG,  GLS,  nor  VMS,  however,  preclude  the  formation  of  over-  and 
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under-shoots  in  the  proximity  of  sharp  gradients  of  the  solution.  For  this 
reason,  discontinuity  capturing  (DC)  techniques,  also  referred  to  as  Spurious 
Oscillations  at  Layers  Diminishing  (SOLD)  methods  are  used  in  combina¬ 
tion  with  SUPG  and  VMS  to  introduce  an  additional  term  to  the  stabilized 
form  of  the  equation.  This  issue  was  treated  for  the  first  time  in  [24],  where 
details  on  how  to  build  the  stabilization  parameter  are  also  given,  and  in  [25] 
for  non-linear  problems.  A  detailed  review  of  most  existing  SOLD  schemes 
can  be  found  in  a  two-part  paper  by  John  and  Knobloch  [26,  27],  where  a 
modification  of  the  discontinuity -capturing  of  Codina  [28]  is  presented  and  is 
shown  to  be  a  promising  option  for  FE  solutions  characterized  by  boundary 
layers. 

All  these  methods  strongly  depend  on  a  parameter  that  will  be  iden¬ 
tified  by  t  throughout  the  paper.  It  will  be  also  referred  to  as  intrinsic 
time.  A  classical  result  for  r  was  obtained  by  Franca,  Frey  and  Hughes  in 
[29]  by  error  analysis.  Their  result  was  reproduced  by  other  authors  using 
different  approaches.  Additional  expressions  for  r  were  found  by  Codina  in 
[30,  31],  by  Codina,  Onate  and  Cervera  in  [32],  by  Harari  and  Hughes  in 
[13],  and  by  Shakib,  Hughes  and  Johan  in  [33],  who  based  the  derivation 
on  the  (discrete)  maximum  principle.  Another  expression  is  due  to  Franca 
and  Valentin  [15]  who  based  their  derivation  on  convergence  and  stability 
analysis.  Starting  with  the  formalization  of  VMS  methods  by  Hughes  [6] ,  r 
has  often  been  derived  using  Green’s  functions,  a  thorough  analysis  of  which 
is  done  by  Hughes  and  Sangalli  in  [34].  Recently,  Houzeaux,  Eguzkitza  and 
Vazquez  [35]  proposed  a  new  way  to  derive  the  approximate  subgrid  scale 
solution,  with  results  that  are  comparable  to  those  of  Hauke  and  Garcia- 
Olivares  in  [36].  In  [37],  Codina  builds  t  using  the  Fourier  analysis  of  the 
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problem;  however,  determining  r  remains  open.  For  this  reason,  we  propose 
t  for  higher-order  spectral  elements  and  use  it  to  construct  an  appropriate 
stabilization  method.  To  further  improve  the  solution,  we  combine  VMS 
and  DC  with  a  method  that  we  here  call  FOS  (for  First-Order  Subcells). 
This  technique  subdivides  a  tensor  product  spectral  element  of  order  p  and 
dimension  d  into  pd  subcells,  and  then  uses  lst-0rder  basis  functions  and 
integration  rules  on  every  subcell  of  the  element. 

1.1.  Main  contribution  of  this  paper 

The  main  problem  that  we  want  to  solve  with  the  work  presented  here 
is  that  of  stabilizing  the  spectral  element  solution  of  the  advection-diffusion 
equation  by  sub-grid  scale  stabilization  techniques  (namely,  VMS),  and  im¬ 
prove  the  solution  by  reducing  the  under-  and  overshoots  that  would  occur 
if  classical  filters  were  to  be  used.  The  definition  and  implementation  of 
t  in  VMS  stabilization  may  greatly  affect  the  result.  We  hence  adapted 
the  method  described  in  [35]  to  compute  r,  and  applied  it  to  spectral  ele¬ 
ments  that  use  Legendre-Gauss-Lobatto  (LGL)  nodes,  and  show  that  this 
technique  can  be  used  in  problems  where  spectral  element  filters  fail.  We 
finally  apply  FOS  (i.e.,  we  lower  the  order  of  interpolation  by  keeping  the 
computational  grid  untouched)  only  where  the  solution  is  characterized  by 
a  propagating  discontinuity.  The  combination  of  VMS,  shock  capturing  and 
FOS  methods  will  be  shown  to  be  an  encouraging  direction  to  take  for  con¬ 
structing  high-order  positive-definite  spectral  element  methods.  To  assess 
the  algorithm,  steady  and  transient  advection-diffusion  problems  are  solved 
on  one-  and  two-dimensional  domains  using  spectral  elements  up  to  order 
16.  We  compare  the  performance  of  the  method  using  VMS  against  those 
obtained  with  previous  classical  spectral  element  schemes. 
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1.2.  Outline  of  this  paper 

The  remainder  of  the  paper  is  organized  as  follows.  After  this  introduc¬ 
tion,  the  numerical  method  and  the  corresponding  stabilization  are  derived 
in  Section  2.  Mass  conservation  properties  are  reported  in  Section  3.  Tests 
to  verify  the  algorithm  and  a  discussion  of  the  results  are  presented  in  Sec¬ 
tions  4  and  5,  respectively. 

2.  Numerical  methods 

Given  the  space  L2  of  real-valued  functions  that  are  square  integrable 
in  a  bounded  domain  12  C  R2  with  boundary  <917,  the  Sobolev  space  Hl  of 
weakly-differentiable  functions  will  be  used.  Specifically,  W  C  H 1  represents 
the  space  of  trial  and  basis  functions  of  the  Galerkin  formulation  to  follow. 
In  L 2  the  inner  product  is  given  by  (•,■),  and  the  2-norm  associated  with 
the  space  is  denoted  by  ||  •  || 2 -  For  simplicity,  we  add  the  property  that 
the  solution  q  vanishes  on  <912;  under  this  assumption,  W  C  Hq.  Given  a 
finite  element  partition  Qh  =  U£i  FQ  of  the  computational  domain  17  into 
nei  high-order  conforming  quadrilaterals  Kj  of  characteristic  length  h,  Wh 
is  the  finite  dimensional  projection  of  W.  The  discrete  weak  form  reduces  to 
the  problem  of  finding  the  function  qh  €  (Wh;  0,  t)  such  that 


(2) 


where  a(-,  •)  is  a  bilinear  form  that  satisfies 


a(fj,q)  =  O,  £(<?))• 


After  integrating  by  parts  and  assuming  homogeneous  Dirichlet  boundary 
conditions,  we  have  that 


Wh,f)=  [  hfdn- 

Jcih 


\h  =  0. 


Remark  1:  If  advection  dominates  diffusion,  unless  h  is  sufficiently 
small  or  the  exact  solution  is  globally  smooth,  the  Galerkin  approxima¬ 


tion  expressed  by  (2)  is  such  that  qh  will  suffer  from  severe/unacceptable 


oscillations  [11,  38].  Furthermore,  if  the  discretization  relies  on  high-order 
methods,  Gibbs  oscillations  may  occur  regardless  of  the  size  of  the  grid.  Dif¬ 
ferent  ways  to  improve  stability  will  be  described  in  the  following  sections. 

2.1.  The  spectral  element  method 

Problem  (2)  is  solved  on  a  grid  of  quadrilateral  elements  of  order  p ,  where 
the  element-wise  solution  qh  is  approximated  by  the  expansion  J2k= l  V;fc(x)  Qk(t) 
on  Np  =  (p  +  l)2  collocation  points  within  the  element.  The  expansion 
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functions  if) are  constructed  as  the  tensor  product  from  the  Lagrange  poly¬ 
nomials  hi(£(x))  and  hj(r](x))  of  order  p  as: 

^fc(x)  =  /ii(£(x))  0  hj(rj(-x)),  y  i,j  =  1,  —,p  +  1-  (4) 

hi(£(x))  and  hj(r](x))  are  the  polynomials  associated  with  the  LGL  points 
and  rjj,  respectively.  The  LGL  points  are  the  zeros  of 

(i-e)p'N(o  =  o 

where  P'N  is  the  derivative  of  the  Nth- order  Legendre  polynomial.  Quadra¬ 
ture  is  performed  on  the  reference  element  Clh  =  [ — 1, 1] 2  with  LGL  points 
that  have  quadrature  weights  u.  Substitution  of  J2k= l  V’fc(x)  int°  the 
weak  form  (2)  yields  the  semi-discrete  (in  space)  matrix  problem 

+  Aqh  +  Bqh  =  0  (5) 

where  q/l  is  the  array  of  the  unknowns  on  the  grid  points,  and  M,  A,  and  D 
are  the  global  mass,  advection,  and  diffusion  matrices,  respectively.  These 
matrices  are  obtained  from  the  direct  stiffness  summation  (DSS)  of  the  ele¬ 
mental  matrices  Mei,  Aei,  and  Dp/  given  by: 


M  ekli=  [  MidVLel 

(6a) 

Jnel 

A %[=  /  u  •  tpi  dttel 

(6b) 

.In*' 

D hi  =  [  vV'tpk  ■  dnet. 

Jrtel 

(6c) 
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By  construction,  the  mass  matrix  M  is  diagonal  (assuming  inexact  integra¬ 
tion)  and  invertible.  This  makes  these  methods  particularly  well  suited  for 
explicit  time  integration  schemes. 

All  the  integrals  defined  above  are  approximated  by  the  quadrature  for¬ 
mula 


r  rl  rl  A  A 

h(')dx=  /  ~  (7) 

J^el  J-lJ-1  1=1  j=l 

where  J  is  the  Jacobian  matrix  associated  with  the  map  between  the  physical 
element  Qh(x,  y )  and  the  reference  element  rj).  The  integration  is  exact 

up  to  polynomials  of  order  2p-l.  The  p+1  LGL  points  lie  along  the  edges 
and  in  the  interior  of  the  elements.  For  more  on  SEM  see,  e.g.,  [39,  40]. 

Integration  in  time  of  (5)  is  performed  with  an  appropriate  strong- 
stability  preserving  (SSP)  time  integrator.  In  particular,  we  use  a  five-stage 
explicit  third-order  Runge-Kutta  method  (RK35)  [41].  SSP  methods  avoid 
the  production  of  additional  oscillations  or  damping. 

From  3rrf-order  and  up,  the  disposition  of  the  nodes  of  spectral  elements 
differs  from  classical  finite  elements  in  the  way  represented  in  Figure  1  for 
a  4ift-order  element. 
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Figure  1:  Nodes  disposition  for  4t,l-order  FE  (left),  and  SE  (right). 


2.2.  Stabilization  techniques 

Aliasing  and  Gibbs  oscillations  render  the  Galerkin  solution  of  (1)  un¬ 
stable  since  unwanted  oscillations  will  pollute  the  numerical  solution.  In 
the  framework  of  spectral  elements,  the  common  strategy  to  control  these 
oscillations  is  the  use  of  anti-aliasing  filters  (See,  e.g.,  [42,  43,  44,  40,  45] 
and  references  therein).  Filtering,  however,  suffers  from  non-positivity  that 
is  unacceptable  in  most  problems  that  involve  transport.  A  suitable  alterna¬ 
tive  to  filters  may  be  the  use  of  a  stabilized  spectral  element  approximation 
based  on  a  residual-based  diffusion-like  term  added  to  the  left  hand-side 
( LHS)  of  (2).  A  stabilization  technique  should  have  the  important  property 
of  consistency  (for  instance,  artificial  diffusion  stabilizes  but  is  not  necessar¬ 
ily  consistent).  Thus,  the  additional  viscous  term  should  vanish  as  the  size 
of  the  element  approaches  zero.  The  stabilized  counterpart  of  (2)  is: 


+  a^h,qh)  +  b^h,qh) 


V^*  €  Wh ,  (8) 


where  b(i/jh,qh)  is  the  stabilization  term.  Different  possible  definitions  of 
b(i/jh,qh )  are  reported  in  Table  1. 
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Table  1:  Stabilized  methods 


Method 

b(i>h,qh) 

cwh) 

AV/HV 

fn,  [V V]  dsr 

£  =  \/aiph 

SU 

fa,  [u  ■  Vq'1]  dfi' 

£  =  u  •  X7i/jh 

SUPG 

fa,  h)r  [dtqh  +  u  ■  Vqh  -  V  ■  {vWqh)  -  f]  dn' 

h> 

II 

e 

<1 

-e- 

3- 

GLS 

fn,  C{^h)T  [dtqh  +  U  ■  Vqh  -  V  •  {vS7qh)  -  /]  dfi' 

£  =  u  ■  Vz/)'1  -  vA4)h 

VMS 

-fa,  [dtqh  +  u  ■  Vqh  -  V  •  (vVqh)  -  f]  dfi' 

£*  =  -u  ■  Viph  -  vAiph 

In  Table  1  and  whenever  they  are  found  below,  the  integrals  on  Q1  are 
defined  by 


/  (-)^'  =  E  /  (O^eZ,  (9) 

Jn>  ^  Jnel 

where  ft'  is  the  union  of  the  element  interiors  only. 

In  the  case  of  AV/HV  and  SU,  b(i/jh,qh)  is  a  function  of  a  constant 
diffusivity  coefficient,  z/,  and  a  =  (3/ 2,  where  j3  is  a  positive  even  power  of 
the  hyper- viscosity  operator  (j3  =  2  yields  the  usual  AV).  Although  HV  is 
scale-selective  (i.e. ,  it  damps  only  higher  frequencies),  it  is  not  consistent, 
nor  is  it  physical.  In  fact,  to  maintain  the  correct  physical  dimensions  of 
the  hyper-viscous  operator,  the  value  of  v  must  be  different  when  different 
a  are  used.  Its  selection  is  hence  not  trivial.  Furthermore,  as  Figure  21  -d 
indicates,  the  diffusion  is  isotropic  and  spatially  homogeneous.  The  operator 
does  not  incorporate  the  problem’s  physics.  On  the  other  hand,  in  the  case  of 
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SUPG,  GLS,  and  VMS,  b(ifh,  qh )  depends  on  what  is  referred  to  as  intrinsic 
time ,  r  (see  Table  1).  r  is  not  constant  or  uniform;  it  depends  on  the 
characteristics  of  the  grid  and  of  the  physics  of  the  problem  (e.g.,  velocity, 
diffusion).  Its  definition  has  major  influence  on  the  accuracy  of  the  solution. 
A  general  method  for  finding  r  does  not  exist  yet.  The  literature  on  the 
optimal  selection  of  r  is  vast,  and  we  refer  to  [26,  27]  for  a  comprehensive 
analysis  of  different  definitions.  Subsection  2.3  and  Appendix  A  show  how 
we  build  t  for  high  order  spectral  elements.  The  method  of  Douglas  and 
Wang  (DW)  [46]  was  omitted  from  Table  1  as  it  can  be  included  into  the 
VMS  method,  although  its  derivation  was  specifically  done  in  the  context 
of  Stokes  problems  rather  than  scalar  advection-diffusion.  SUPG  and  VMS 
will  be  described  in  the  following  subsections. 

2.2.1.  Streamline-upwind/ Petrov-  Galerkin  (SUPG) 

The  SUPG  method  was  designed  by  Brooks  and  Hughes  [9]  and  was 
later  generalized  for  multidimensional  problems  by  Hughes  and  Mallet  [47]. 
It  is  a  consistent  alternative  to  the  artificial  diffusion  approach  or  to  the 
overly  diffusive  streamline  upwind  (SU)  method.  Its  use  has  been  ubiquitous 
in  the  solution  of  transport  problems  by  the  finite  element  method  (See, 
e.g.,  [48,  29,  49,  50,  51]).  The  application  of  this  strategy  to  higher-order 
schemes  was  first  tested  for  spectral  methods  by  Canuto  and  coworkers  in 
[17,  18,  19,  52],  and  later  by  Hughes  and  coworkers  in  [21]  using  non-uniform 
rational  B-splines  (NURBS).  In  this  paper  we  show  its  properties  when  used 
with  high-order  spectral  elements.  SUPG  is  a  Petrov-Galerkin  method  in 
that  it  does  not  assume  that  the  basis  and  test  functions  live  in  the  same 
space.  We  introduce  the  additional  space  of  test  functions  wh  defined  by 
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S <h  =  [wh  :  wh  =  ifh  +  ru  •  Viph  :  V> h  €  . 


We  have  the  problem  of  finding  the  function  €  (TV^-;  0,  t)  such  that 


+  tu  •  V^, 


cV\ 
"dt"  ) 


+a('iph+Tu-Viph,qh) 


('iph+Tu-\7iph,  f) 


Vifh  €  Wh. 

(10) 


Rearrangement  of  (10)  yields 


+  a(ifh,qh)-(iph,f)+b(ifh,qh)  =  0 

■  v - ^ 

Galerkin 


€  Wh  (11) 


where 


b{'iph,  qh) 


dqh 

~dt 


+  u  •  X7qh 


V  •  (vVqh)  -  f 


t  u  •  vv>h  dnh 


(12) 


is  the  stabilizing  term.  In  (12),  dtqh  +  u-  S7qh  —  V  •  (vS7qh)  —  /  is  the  residual 
of  the  governing  equation,  and  r  is  the  stabilization  parameter  whose  con¬ 
struction  is  reported  below  in  the  specific  context  of  variational  multiscale 
stabilization. 

2.2.2.  The  variational  subgrid  scale  formulation  (VMS) 

Let  Wh  be  the  space  of  resolved  scales  and  let  W  be  a  space  that  com¬ 
pletes  Wh  in  W  and  that  we  will  call  the  space  of  subgrid  scales.  VMS 
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relies  on  the  decomposition  W  =  Wh  ©  W,  where  ©  is  the  overlapping  sum 
decomposition  of  the  two  spaces.  The  elements  that  belong  to  W  are  q  and 
ip  and  are  such  that  q  =  qh  +  q  and  ip  =  iph  +  ip.  Using  such  decomposition 
of  q  and  ip  and  anticipating  that  we  will  consider  q  to  be  quasi-static  [37], 
we  re-write  (2)  as 

(?h  +  $’l(r)  +  a(^  +  i<^  +  g)  =  (V’'l  +  ^/)  €  wh,ip  €  w. 

(13) 

By  virtue  of  the  linear  independence  of  iph  and  ip  we  can  first  take  ip  =  0 
and  then  iph  =  0  and  find  the  split  problem: 


(^h,^j+a(iPh,qh)  +  a(iPh,q)  =  (iPh,f)  ViPh  €  Wh  (14a) 

($,  j  +  «(V’>  Qh)  +  <*($,  Q)  =  V-0  e  W.  (14b) 

We  make  the  assumptions  that  ip(d£l)  =  0  and  q(d£l)  =  0,  and  that 
ip(dQei)  =  0. 

Following  [6],  in  (14)  we  integrate  by  parts  the  bilinear  forms  that  depend 
on  the  subgrid-scales  and  write  the  following: 


a(iph,q)  =  (C*iph,q) 
a(iP,qh)  =  (iP,Cqh ) 
a(4>,q)  =  {ip,Cq) 


Viph  €  Wh  and  q  €  W 
V-i/S  €  W  and  qh  €  Wh 
yip  €  W  and  q  €  W, 
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where  C*  is  the  dual  (or  adjoint)  of  C.  The  following  system  is  then  found: 


Uh,^+a^h,qh)  +  (/:*^h,q)  =  ^hJ)  €  Wh  (15a) 

($,  ^-)  +  +  $,£q)  =  (i’j)  V'ipeW,  (15b) 

where 

(C*^h,q)=  [  £*(il>h)qdn', 

.nv 

('4>Xqh)=[  i>C(qh)d£i', 

■nv 

(ti,Cq)  =  j  %i)C(q)dQ! . 

.nv 


Observations  on  time-dependent  subgrid-scales:  The  time-dependent 
approximation  (13)  would  include  a  contribution  from  the  time  evolution  of 
the  subscales  given  by  dtq  if  the  the  hypothesis  of  quasi-static  subscales  (i.e. 
dtq  ~  0)  had  not  be  considered.  Under  this  hypothesis,  the  contribution 
from  the  subgrid  scales  only  appears  in  the  steady  part  of  the  Galerkin  ap¬ 
proximation.  If  a  sufficiently  small  time-step  is  used  with  an  explicit  time 
integrator,  we  do  not  lose  accuracy  with  the  quasi-static  hypothesis.  With 
the  use  of  large  time-steps  with  semi- implicit  time  integrators  in  atmospheric 
simulations,  tracking  of  the  subscales  is  hence  needed.  This  issue  is  reserved 
for  future  work  by  the  authors. 
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Remark  2:  It  must  be  pointed  out  that,  for  the  solution  of  the  scalar 
transport  equation,  the  expressions  for  SUPG  and  VMS  are  the  same. 

Approximation  of  the  sub-grid  scales.  The  unresolved  quantity  q  has  not 
been  defined  in  any  way  yet.  In  the  following,  the  subscales  are  constructed 
by  algebraic  approximation  following  the  technique  described  by  [35]  for  lin¬ 
ear,  quadratic,  and  cubic  elements.  They  take  W  as  the  space  of  vanishing 
functions  on  the  boundaries  of  each  element.  These  are  called  bubble  func¬ 
tions  (see  [53,  54]).  By  incorporating  the  time-dependent  term  of  Eq.  (15b) 
within  the  second  inner  product  and  by  re-arranging  the  terms,  the  equation 
for  the  subgrid  scales  q, 


ipC(q)  <m' 


[£(<?)  -  /]  <m'  g  w, 


(17) 


is  found.  The  equivalent  strong  form  of  (17)  is 


£(?) 


C-{qh)  ~  f 


~R(qh), 


(18) 


where  R(qh )  is  the  residual  of  the  equation  for  the  resolved  scale.  For  the 
solution  of  (18)  we  first  define  q  as  a  function  of  the  bubbles  6(x);  we  have 
q  =  — b{x)R{qh )  that  is  plugged  into  (18)  to  find  the  differential  problem 

£(-b(x)R(q*))  =  -R(qh).  (19) 


Eq.  (19)  is  solved  for  b  with  Dirichlet  boundary  conditions  b(x i)  =  0  and 
b(x 2)  =  0  on  every  element  of  length  h  =  \x\  —  X2\-  By  thinking  that  R(qh)  is 
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always  known  from  the  previous  time-step,  we  can  consider  it  as  a  constant 
and  can  hence  be  taken  out  of  the  £(•)  operator  so  that  C(b(x))  =  1.  For 
the  one-dinrensional  steady-state  advection-diffusion  equation1  the  problem 
is 


C(b(x ))  =  ubx(x )  -  v  bxx(x)  =  1, 


(20) 


with  exact  solution 


_  x  hi-  exu!v 

U  +  U  ehu/v 


(21) 


Now  that  we  computed  the  bubble  functions  along  the  element  lenght, 
we  construct  the  stabilization  parameter  r  as  the  mean  value  of  b(x)  along 
each  element.  We  have  that 


T  = 


(22) 


Integration  of  b(x)  (21)  in  the  interval  [x \ ,  x'2]  =  [0,  h\  yields  the  following 

lrrhe  bubbles  are  computed  only  once,  out  of  the  time  loop;  only  if  they  depend  on  non¬ 
constant  coefficients,  they  are  computed  once  per  time-step.  In  either  case,  the  problem 
to  be  solved  is  not  time-dependent. 


19 


definition  of  r: 


T  =  Yu  (coth(Pei) '  •  <23> 


where 


Pek 


uh 

~2u 


is  the  element  Peclet  number. 

We  have  derived  all  the  ingredients  to  define  the  sub-grid  scales  q  as  the 
algebraic  approximation  [6] 


q  =  -rR{qh).  (24) 

Eq.  (15b)  was  used  as  the  starting  point  to  approximate  q.  Now,  by 
plugging  (24)  into  (15a),  the  VMS  stabilized  Galerkin  method  is  found  and 
expressed  as  follows: 

Find  qh  €  Wh  such  that 


(^h,^pj+a^h,qh)-(f^h)- J^jr*^h)TR(qh)dQ'  =  0  €  Wh. 

(25) 

Eq.  (25)  differs  from  Eq.  (2)  by  the  additional  term  that  models  the  subgrid 
scales.  The  extra  term  is  the  viscous-like  contribution  that  stabilizes  the 
equation. 
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2.3.  t  for  spectral  elements 

For  linear  elements,  h  =  \x\  —  X2I  is  simply  the  length  of  the  element. 
For  higher-order  finite  elements,  h  becomes  a  fraction  of  the  total  element 
size  if  the  internal  nodes  are  equi-spaced.  In  the  case  of  spectral  elements, 
where  the  LGL  points  are  unevenly  distributed,  the  integral  is  computed  by 
using  h  as  the  local  distance  between  two  consecutive  points.  The  stabiliza¬ 
tion  parameter  r  is  built  inside  the  element  as  a  function  of  the  bubbles  on 
every  segment  delimited  by  two  consecutive  nodes.  This  means  that  equa¬ 
tion  C(b(x))  =  1  is  solved  on  every  sub-element  by  applying  homogenoeus 
Dirichlet  b.c.  at  the  sub-element  boundaries.  For  example,  for  a  second- 
order  element  with  one  internal  node  we  would  solve  C(b(x))  =  1  on  the  two 
segments  [x\,x?\  and  [3:2, 3:3],  respectively,  by  applying  homogeneous  b.c.  as 
b(x  1)  =  0,  6(3:2)  =  0  and  6(3:2)  =  0,  6(3:3)  =  0.  With  this,  two  r’s  would  be 
computed  as 


T„*+1  = 


Xlgl(i  +  1)  -  Xlgl{i)  J xlgl (2) 


b(x)  dx, 


(26) 


where  xigi(i  +  1)  and  xigi(i )  are  the  coordinates  of  two  consecutive  LGL 
points. 

The  most  simple  (but  not  unique)  way  of  proceeding  is  that  of  taking 
the  average  value  of  all  the  sub-r’s  as  the  value  of  the  full  element  r.  This 
method  is  used  in  this  paper. 

The  uneven  spacing  of  the  element  nodes  is  the  major  difference  with 
respect  to  the  definitions  derived  in  previous  studies.  In  this  case  the  intrin¬ 
sic  time  is  non-uniform  along  the  element,  as  it  appears  in  Figure  2,  where 
the  bubbles  and  corresponding  r’s  are  displayed  for  an  element  of  order 
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Figure  2:  Bubbles  b(x)  and  r  for  a  7th-order  unitary  spectral  element.  The  biggest  bubble 
in  the  plot  is  the  bubble  that  a  linear  element  would  have. 

7.  In  Appendix  A  we  report  the  explicit  expression  to  compute  along  a 
high-order  spectral  element. 

2-4-  Spurious  oscillations  at  layers  diminishing  (SOLD)  methods 

Methods  in  the  form  of  (25)  may  produce  overshoots  and  undershoots 
in  the  proximity  of  an  internal  or  boundary  layer.  These  unwanted  oscilla¬ 
tions  can  be  suppressed,  without  affecting  the  global  solution,  by  adding  an 
additional  diffusive  term  of  the  form 

Wfc,f  Vqh),  (27) 


where  consistency  must  be  respected  through  a  proper  construction  of  f. 
We  would  like  to  have  a  method  that  does  not  modify  the  diffusion  in  the 
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streamline  direction  since  that  is  already  accounted  for  by  the  stabilization 
term,  but  also  avoids  overdamping  in  the  crosswind  direction.  The  compre¬ 
hensive  set  of  tests  performed  by  John  and  Knobloch  reveals  that  Codina’s 
[28]  is  among  the  best  methods  that  satisfy  these  conditions  when  used  with 
finite  elements.  In  [28],  f  is  defined  by: 


1 

f  =  -max 
2 


2v 

U|||6'fc 


hk 


\R(qh)\ 

liv^n 


u  0  u 


u 


(28) 


where  C is  a  constant,  U||  is  the  velocity  component  in  the  streamline  direc¬ 
tion,  and  0  indicates  a  tensor  product.  Codina  suggests  C  =  0.7  for  linear 
and  bilinear  elements,  and  C  =  0.35  for  quadratic  and  biquadratic  elements. 
However,  for  higher  order  elements  using  LGL  points,  we  found  that  the  best 
results  were  obtained  by  setting  C  =  1,  as  long  as  hk  is  selected  properly  in 
the  construction  of  both  f  and  Tfc. 


An  alternative  to  (27)  comes  from  Johnson,  Schatz,  and  Wahlbin  [55] 
who  defined  the  following: 

(fux  •  VY^,  ux  •  Vqh)  ,  ux  =  (29) 

In  the  current  work,  (29)  gives  better  results  than  (27),  and  was  then  used 
throughout.  The  results  obtained  with  this  technique  are  labeled  with  DC 
for  Discontinuity  Capturing. 

2.5.  First-Order  Subcells  (FOS) 

FOS  is  one  additional  tool  that  can  further  help  the  suppression  of  Gibbs 
oscillations.  The  concept  is  simple  and  is  easily  coded  on  structured  grids.  If 
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a  solution  has  large  gradients,  the  algorithm  needs  to  identify  the  elements 
where  the  large  gradients  occur,  and  project  the  solution  scheme  to  a  1st- 
order  space.  The  gradient  is  sought  with  a  proper  error  estimator.  The 
simple  physics  of  the  advection-diffusion  problems  discussed  below  allows 
for  the  energy-norm  of  the  gradient  of  the  solution  to  be  a  sufficiently  good 
estimator  for  the  current  study.  As  it  is  defined  in  this  study,  the  error 
estimator  depends  on  a  parameter,  e,  that  may  be  a  function  of  the  numerical 
settings  (e.g.,  grid  resolution,  time  step).  This  point  must  be  considered  in 
the  construction  of  FOS  and  in  the  selection  of  the  error  estimator.  We  did 
not  explore  this  further  in  this  study,  although  it  is  a  very  important  issue 
for  the  best  performance  of  FOS. 

Algorithm  1  is  a  simple  implementation  of  this  concept  within  our  code. 
We  present  the  pseudo-code  below  for  the  sake  of  clarity.  The  method  was 
applied  to  a  two-dimensional  advection-diffusion  problem  with  internal  and 
boundary  layers  in  a  skew  velocity  field.  Results  are  shown  in  Figures  9-12 
and  14-18.  A  detail  of  Figures  18(a,6)  is  presented  in  Figure  19. 

In  the  tests  that  use  Algorithm  1,  e  was  set  to  0.5. 

3.  Mass  conservation 

For  problems  in  geophysical  fluid  dynamics  and,  more  specifically,  in  at¬ 
mospheric  simulations,  mass  conservation  of  tracers  is  an  important  ingre¬ 
dient.  In  this  section  we  address  this  issue  and  illustrate  how  the  algorithm 
reported  in  this  paper  behaves  in  this  respect.  Under  the  suitable  hypothesis 
of  a  divergence-free  flow,  Eq.  (1)  for  the  transport  of  the  mixing  ratio 
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for  iel  =  1  to  nelem  do 

//  Check  if  the  element  contains  a  sufficiently  large  gradient: 
if  iel  is  s.t.  || \7qh  II2  >  e  then 

//  Treat  element  iel  as  a  sub-domain  made  of  ( ngl  —  1)  x  ( ngl  —  1) 
sub-elements  (isubel) : 
for  isubel  =  1  to  pd  do 

Create  mass  and  rhs  using  l,si-order  basis  functions  and 
integration  rule 

end  for 

else 

Create  rhs  for  the  high-order  spectral  element. 

end  if 

end  for 

Algorithm  1:  Compute  the  lst-order  rhs 
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where  ptracer  and  p  are,  respectively,  the  densities  of  the  tracer  and  of  the 
advecting  fluid,  is  derived  from  the  equation  of  conservation  of  mass  of  the 
tracer2 


—  +  V  •  (u pq)  =  0  (30) 

by  elimination  of  p  from  the  conservation  law  (continuity  equation) 

+  V  ■  (up)  =  0.  (31) 


In  (30)  and  (31)  the  quantities  that  are  conserved  are  ptracer  =  PQ  and  p, 
but  not  q.  In  the  case  of  p(t)  =  p(t  =  0)  =  constant  in  a  non-divergent  flow, 
Eq.  (1)  is  equivalent  to 


+  V  •  (uq)  =  0.  (32) 

This  allows  the  equal  treatment  of  pq  and  q  in  the  modeling  of  the  physi¬ 
cal  system  at  hand  [56].  Regardless  of  the  poor  conservation  properties  of 
transport  expressed  in  advective  form,  Eq.  (1)  is  of  common  use  within  at¬ 
mospheric  models  (e.g.,  [1,  56,  57]).  In  this  paper  we  track  mass  loss  during 
the  simulations  to  evaluate  the  amount  of  mass  loss  that  we  would  run  into 
if  using  SEM+VMS  knowing  that  no  method  will  conserve  for  this  equation 
since  it  is  not  in  conservation  form  or  a  conservation  law. 

2 For  simplicity,  the  diffusion  and  source  terms  were  dropped. 
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For  this  computation,  we  use  the  advection  of  a  square  wave  in  a  periodic 
channel  and  report  the  results  in  paragraph  Tr2-2D  below. 

4.  Numerical  testing 

The  algorithms  discussed  throughout  this  paper  are  tested  by  using  stan¬ 
dard  one-  and  two-dimensional  problems.  The  problems  are  organized  ac¬ 
cording  to  the  nomenclature  listed  below: 

•  ID  Steady-state  homogeneous  advection-diffusion  (St- ID) 

•  ID  Steady-state  advection-diffusion  with  source  ( St-ID-S) 

•  2D  Steady-state  advection-diffusion  with  internal  and  boundary  layers 
(St- 2D). 

•  2D  Time-dependent  advection-diffusion  with  “L”-shaped  discontinuity 

(: Trl-2D ). 

•  2D  Time-dependent  advection  of  a  sharp  tracer  in  a  doubly  periodic 
channel  (Tr2-2D). 

•  2D  Smooth  solid-body  rotation  -  Convergence  study  (2D  Smooth  solid- 
body  rotation ). 

•  Pseudo-3D  advection  in  a  neutrally  stratified  atmospheric  flow  (Atmo- 
3D). 

St-ID.  One-dimensional  steady-state  advection-diffusion.  The  tracer  qh  is 
propagated  with  constant  velocity  u  =  1  m  s_1  and  diffusivity  u  =  1/512  m2  s-1 
first  on  two  elements  of  order  p  =  10  (Figure  3),  and  then  on  four  elements 
of  order  p  =  12  (Figure  4).  The  domain  is  the  line  segment  D  =  [—1,1]  with 
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Dirichlet  boundary  conditions  qh(— 1)  =  0  and  qh(l)  =  1.  We  compared 
the  filtered  (top  row)  against  the  stabilized  solution  (bottom  row)  and  ob¬ 
serve  a  decrease  of  oscillations  and  undershoots.  Also,  at  higher  order  and 
finer  resolution,  the  capabilities  of  the  filter  are  clearly  being  challenged 
by  the  presence  of  the  boundary  layer  at  x  =  1.  At  the  same  time,  small 
oscillations  near  the  nodes  of  the  element  by  the  boundary  layer  are  not 


completely  suppressed  by  the  stabilized 
localized  smoothing  is  sought. 


X 


(a)  Filter.  2  el,  p  =  10 


x 


(b)  VMS.  2  el,  p  =  10 

Figure  3:  St-ID :  v  =  l/512m2s_1.  The  ex¬ 
act  solution  is  dashed.  The  circles  indicate 
the  grid  points. 

28 


method  either;  hence,  additional 


X 


(a)  Filter.  4  el,  p  =  12 


x 


(b)  VMS.  4  el,  p  =  12 

Figure  4:  St-ID:  v  =  1/512  m2  s_1.  The  ex¬ 
act  solution  is  dashed.  The  circles  indicate 
the  grid  points. 


St-ID-S.  Steady  state  advection-diffusion  with  source  term  /  =  1.  qh  is 
propagated  with  constant  velocity  u  =  lm  s-1  and  two  different  diffusivities: 
v  =  5  x  10-3  m 2  and  v  =  5  x  10~2  m2  s  .  The  domain  is  the  line  segment 
17  =  [0, 1]  and  homogeneous  Dirichlet  boundary  conditions  are  imposed. 
The  domain  is  subdivided  into  two  elements  of  order  p  =  16  and  runs  are 
compared  using  filtered  SE  (top  row  in  Figures  5  and  6),  and  VMS  (bottom 
row).  We  observe  a  very  similar  behavior  of  the  solution  among  the  two 
different  cases  in  the  smooth  problem  (Figure  5).  The  results  are  comparable 
to  the  ones  obtained  by  Houzeaux  et  al.  with  their  r  for  quadratic  and  cubic 
elements  in  [35]. 
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(a)  Filter  (a)  Filter 
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X 

(b)  VMS 

(b)  VMS 

Figure  5:  St-ID-S:  v  —  5  x  10 ~3m2s_1.  2 
16th-order  elements.  The  exact  solution  is 
dashed.  The  circles  indicate  the  grid  points. 


Figure  6:  St-ID-S:  v  =  5  x  10-2  m2s_1.  2 
16th-order  elements.  The  exact  solution  is 
dashed.  The  circles  indicate  the  grid  points. 


St- 2D.  Standard  steady  advection-diffusion  skewed  to  the  mesh  (e.g.,  [28]): 
a  discontinuity  is  propagated  with  constant  velocity  u  =  (1,-2 )ms~1  and 
diffusivity  v  =  10~8  m2  s-1  in  the  unit  square  12  =  [0, 1]  x  [0, 1].  The  initial 
configuration  is  shown  in  Figure  7. 
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qh  =  1 


B 


%  =1 


(1,-2) 


qh  =o 


qh  =° 

Figure  7:  St-2D :  initial  configuration  of  the  steady-state  problem 


Dirichlet  boundary  conditions  are  prescribed: 


1 


if  y  =  1, 


q 


h 


1  if  x  =  0  and  y  >  0.7, 


0  otherwise. 


In  Figures  8-12  we  illustrate  the  run  of  the  same  case  with  different  number 
of  elements  and  order  of  the  interpolating  polynomials.  For  direct  compari¬ 
son  of  our  solution  with  the  ones  in  the  existing  literature  of  finite  elements, 
we  first  run  the  test  with  linear  elements  (p  =  1),  and  present  the  results  in 
Figure  8.  The  multiscale  solution  of  this  problem  (see  Figure  8  a)  shows  im¬ 
portant  boundary  and  internal  layers  that  are  damped  with  the  discontinu¬ 
ity  capturing  techniques  of  Section  2.4.  The  application  of  the  discontinuity 
capturing  (DC)  scheme  greatly  improves  the  solution  and  yields  monotonic¬ 
ity  (see  Figure  86).  In  Figure  9  we  maintained  the  same  number  of  nodes 
of  the  previous  run,  but  increased  to  4th  the  order  of  interpolation  to  as¬ 
sess  the  algorithm  in  the  context  of  this  paper  (i.e.  50  elements  of  order 


31 


4  were  used  instead  of  200  elements  of  order  1).  The  similar  behavior  of 
the  solution  with  respect  to  the  lst-order  polynomial  run  suggests  that  the 
residual-based  methods  as  implemented  in  this  study  may  not  be  sensitive  to 
the  distribution  of  the  interpolation  nodes  within  the  elements  edges.  As  it 
appears  in  Figures  9  c,  the  behavior  is  completely  analogous  to  the  previous 
run.  However,  monotonicity  is  lost  in  two  singular  nodes:  with  reference  to 
Figure  9  c,  the  4<fe-order  solution  is  smooth  and  monotone  everywhere  except 
for  the  nodes  represented  as  points  A  and  B  in  Figure  7.  This  is  not  surpris¬ 
ing:  at  A  and  B  the  tracer  is  leaving  the  boundary  with  a  skew  angle;  an 
incorrect  imposition  of  boundary  conditions  at  these  nodes  may  be  causing 
the  problem.  The  numerical  singularity  at  this  points  should  be  addressed 
but  it  will  not  be  done  in  the  current  work.  These  are  fully  suppressed  by 
applying  the  FOS  algorithm  described  above,  as  it  is  shown  in  Figure  9 d. 

Decreasing  the  number  of  computational  nodes  by  doubling  the  order 
from  4  to  8  and  setting  the  number  of  elements  to  10  in  x  and  z,  even  with  a 
discontinuity  capturing  term,  the  solution  starts  to  lose  nronotonicity.  This 
appears  in  Figures  11  and  12,  where  extrema  get  larger  than  in  the  previous 
cases.  This  problem  shows  that  the  construction  of  the  stabilizing  parameter 
r  should  include  information  on  the  order  of  the  interpolating  polynomial. 

For  a  better  view  of  the  problem,  in  Figures  10  and  12  we  present  a  ver¬ 
tical  slice  of  the  solution.  The  boundary  layers  are  evident.  Their  damping, 
however,  is  clear  if  VMS,  DC,  and  FOS  are  applied. 
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(a)  VMS  (b)  VMS  +  DC 

Figure  8:  St- 2D:  steady-state  solution  on  200x200  lst-order  elements.  (For  plotting  only, 
the  data  are  interpolated  to  a  50  x  50  node  grid  using  Octave  [58]). 
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(Mm,  Max)  =  (-5,509e-3, 1.026)  (Min,  Max)  =  (-2.161e-07,  1) 


(c)  VMS  +  DC  (d)  VMS  +  DC  +  FOS 

Figure  9:  St-2D :  steady-state  solution  on  50x50  4th-order  elements. 
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Figure  10:  St-2D\  steady-state  solution  on  50x50  4tft,-order  elements.  Vertical  slice  at 
z  =  0.3 
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(a)  Filter 


(b)  VMS 


3.5 


Figure  12:  St-2D\  steady-state  solution  on  10x10  8tft,-order  elements.  Vertical  slice  at 
z  =  0.3 


Trl-2D.  Transient  advection-diffusion  of  an  L-shaped  discontinuity  in  a  flow 
where  v  =  1CP6  m2  and  the  velocity  u  of  magnitude  |u|  =  0.5\/2  ms-1 
is  at  45°  with  respect  to  the  axis  ( x,z ).  The  initial  configuration  is  shown 
in  Figure  13. 
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Oh  =  i  b  qh  =  o 


Figure  13:  Trl-2D:  initial  configuration  of  the  L-shaped  problem 

The  convex  shape  of  the  sharp  discontinuity  makes  this  problem  more 
challenging  than  the  previous  case  [59],  and  is  chosen  to  analyze  robustness 
and  accuracy  of  the  algorithm.  Runs  were  performed  at  two  different  resolu¬ 
tions  and  two  different  orders  of  interpolating  polynomials.  In  particular  we 
have:  approx.  100  points  per  side  using  25x25  4th  —  order  elements  (Figure 
14),  and  12x12  8tft-order  elements  (Figure  15);  and  approx.  200  points  per 
side  using  50x50  4<ft-order  (Figure  17),  and  25x25  8th-order  (Figure  17).  In 
the  figures,  Filter  means  that  the  SEM  solution  was  filtered  at  every  time- 
step.  VMS  and/or  DC  indicate  that  the  SEM  solution  is  stabilized  by  the 
VMS  with  or  without  a  discontinuity  capturing  term  (DC).  VMS  +  DC  + 
FOS  indicates  the  contribution  of  FOS  as  well.  Positivity  is  not  preserved 
in  the  solution  obtained  with  a  filter.  The  sharp  front,  in  fact,  makes  the  fil¬ 
ter  inappropriate.  However,  similarly  to  the  steady  advection-diffusion  test 
St-2D,  the  VMS-stabilized  solution  of  this  problem  is  characterized  as  well 
by  the  formation  of  internal  layers  that  run  along  the  edges  of  the  tracer  in 
the  direction  of  the  flow  (See,  e.g.,  Figure  14b),  and  VMS  is  not  sufficient  to 
preserve  monotonicity  unless  it  is  supplemented  by  the  additional  DC  term 
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defined  in  (29).  This  effect  is  displayed  in  Figures  14,15,16,  and  17. 

The  consideration  made  for  problem  St-2D  on  the  singular  peaks  that 
form  at  the  nodes  where  the  tracer  leaves  the  boundary  at  an  angle,  applies 
here  at  nodes  A  and  B  of  Figure  13.  This  is  visible  in  Figure  18  obtained 
by  slicing  the  tracer  along  z  =  0  in  Figures  16  and  17,  respectively.  The 
problem  is  solved  by  the  application  of  FOS. 

As  the  order  of  interpolation  is  increased  from  4th  to  8th,  the  smooth 
solution  begins  to  lose  positivity.  As  interpreted  for  St- 2D,  the  solution  is 
clearly  being  affected  when  the  interpolation  nodes  are  densely  clustered 
towards  the  boundaries  of  the  elements,  as  is  the  case  for  higher  order. 
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(a)  Filter 


Figure  16:  Trl-2D :  4' 


order  50x50.  t 


■HUH 


(a)  Filter 


order  25x25.  t 


(a)  4th  —  order  50x50 


(b)  8th  -  order  25x25 
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Figure  18:  Trl-2D\  Vertical  slice  at  z  =  0.0. 


X 


(a)  4th  —  order  50x50 


(b)  8th  -order  25x25 
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Figure  19:  Trl-2D:  detail  of  Figure  18.  Region  with  undershoots.  Vertical  slice  at  2  =  0.0. 


Tr2-2D.  Linear  advection  of  a  2D  square  wave  along  x  in  the  periodic 
domain  D  =  [0,1]  x  [0,1]:  the  tracer  is  transported  with  velocity  u  = 
(1/2,  0)  ms-1  for  one  periodic  revolution  along  x.  The  initial  concentration 
qh  =  1  is  centered  at  (xc,zc)  =  (0.5,  0.5)  (Figure  20).  The  computational 
finite  domain  consists  of  11  x  11  quadrilaterals  of  order  11. 


Figure  20:  Tr2-2D\  initial  configuration  of  the  pure  advection  problem. 

As  in  the  steady  case,  Figures  21-23  display  improvement  of  the  solution 
in  terms  of  nronotonicity  when  the  VMS  method  is  used  instead  of  the 
filter.  The  combination  of  VMS  and  filtering  is  not  recommended  (result 
not  shown);  although  VMS  alone  controls  the  over-  and  under-shootings 
along  the  streamlines,  the  addition  of  the  filter  at  the  end  of  every  time  step 
degrades  positivity  in  the  neighborhood  of  large  gradients. 

In  Figures  22  and  23  we  present  the  streamline  and  crosswind  sections 
of  the  solution  obtained  by  slicing  the  tracer  along  z  =  0.5  and  x  =  0.5, 
respectively.  Unlike  the  previous  problems  characterized  by  internal  and 
boundary  layers,  for  pure  advection  the  VMS  preserves  the  maximum  and 
minimum  concentrations  q^ax  =  1,  and  q^in  =  0  and  is  free  of  spurious 
oscillations.  As  a  point  of  comparison,  we  present  the  result  of  classical 
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artificial-viscosity  in  Figures  21-23 d. 


(Min,  Max)  =  (-0.1 191,  1.11) 


(a)  Galerkin 


(Min.  Max)  =  (-4.125e-06, 1) 


(Min,  Max)  =  (-0.2023,  1 .097) 

(b)  Filter 


(Min.Max)  =  (0.0,  0.9482) 
(d)  AV/HV  V  =  0.001  m2  s"1 


(c)  VMS 

Figure  21:  Tr2-2D:  Surface  plot  of  the  concentration  field:  At  =  0.001s  (except  for  HV: 
At  =  0.0002s),  11x11  elements  with  11th  order  polynomials.  Results  at  t  =  2.0  s  (after 
1  periodic  revolution  along  x). 
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(a)  Galerkin  (b)  Filter 


(c)  VMS  (d)  AV/HV  v  =  0.001  m2  s"1 


Figure  22:  Tr2-2D:  Streamline  cut  at  0.5  m  in  the  y-direction.  At  =  0.001s,  11x11 
elements  with  11th  order  polynomials.  Results  at  t  =  2  s  (after  1  periodic  revolution 
along  x).  Solid  line  indicates  the  computed  solution.  The  dashed  line  is  the  analytic 
solution. 
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(c)  VMS 


(d)  AV/HV  v  =  0.001  m2  s"1 


Figure  23:  Tr2-2D:  Crosswind  cut  at  0.5  m  in  the  x-direction.  At  =  0.001s,  11x11 
elements  with  11th  order  polynomials.  Results  at  t  =  2  s  (after  1  periodic  revolution 
along  x).  Solid  line  indicates  the  computed  solution.  The  dashed  line  is  the  analytic 
solution. 


Testing  mass  conservation.  Because  of  the  periodic  boundary  conditions 
applied  here,  we  compute  mass  conservation  properties  for  this  test.  At 
every  time-step,  the  total  mass  loss  of  pq  (for  p  =  1)  is  computed  as 


Mlossit) 


fn(pq(t)  -  M* q))  dn 

In  PQ(t o)  dtt 


(33) 


where  17  is  the  domain  volume  and  to  indicates  the  values  at  the  initial 
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time.  Figure  24  shows  the  evolution  of  the  mass  loss  that  occurs  during 
100  revolutions  around  the  periodic  channel  of  Test  Tr2-2D.  100  revolutions 
happen  in  100  s. 


Figure  24:  Tr2-2D\  evolution  of  the  mass  loss  during  100  s,  or  100  revolutions  of  the  square 
wave  around  the  periodic  channel  of  Figure  20. 

Although  in  Figure  24  there  seems  to  be  an  asymptotic  trend  to  an  up¬ 
per  bound,  this  is  obtained  at  the  expenses  of  accuracy  during  long  runs. 
Regardless  of  the  type  of  equations  (conservative  or  non-conservative),  the 
method  here  proposed  is  certainly  unable  to  retain  all  mass.  Because  of 
this,  at  this  stage  we  can  only  think  of  applying  this  technique  to  short  term 
weather  forecast  but  not  climate.  This  is  the  first  application  of  VMS  and 
DC  to  Spectral  Elements  to  solve  the  advection  equation;  in  the  future  we 
will  work  on  a  fix  to  this  problem  for  better  (or  total)  mass  conservation.  A 
first  improvement  of  accuracy  for  long-time  runs  may  be  achieved  using  or- 
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thogonal  sub-grid  scales  (OSS)  as  proposed  by  [31].  Further  analysis  should 
be  done. 

2D  Smooth  solid-body  rotation.  The  smooth  solid-body  rotation  test  with  a 
smooth  function  is  used  for  a  grid  convergence  study  [60,  61].  A  Gaussian 
hill  qh  =  exp  [—5  ((x  —  xc )2  +  (y  —  yc )2)]  is  originally  centered  in  (xc,  yc)  = 
(0, 0)  in  a  periodic  domain  17  =  [— n,  7r]  x  [— 7r,  7t]  with  prescribed  velocity 
(■ u,w )  =  (— Try,  nx).  Convergence  is  computed  after  one  full  revolution 
(t  =  2  s).  The  normalized  standard  L2  error  is  computed  with  respect  to 
the  exact  solution  qe  =  qh(x,y,t  =  0)  using  Nei  €  { 102 , 202, 402}  and 
polynomials  of  degree  4  and  8.  Figure  25  shows  the  h— error.  The  original 
data  are  reported  in  Table  2. 


Figure  25:  2D  Smooth  solid-body  rotation :  Log-log  plots  of  the  normalized  L2  error  vs. 
Nei  using  VMS  or  a  Filter. 


49 


Table  2:  2D  Smooth  solid-body  rotation :  normalized  1/2  error  vs.  Nei.  Convergence  rate 
of  every  setting  is  reported  on  the  last  row  of  the  table. 


Net 

VMS  4th 

FILTER  4th 

VMS  8th 

FILTER  8th 

100 

0.7083E-02 

0.1793E+00 

0.3175E-05 

0.2569E-03 

400 

0.3280E-03 

0.2422E-01 

0.2305E-07 

0.4927E-05 

1600 

0.1982E-04 

0.4127E-02 

0.5473E-10 

0.9388E-07 

Convergence  rate: 

4.2406 

2.7206 

7.9120 

5.7091 

The  experiment  indicates  that  VMS  does  not  affect  the  rate  of  con¬ 
vergence  of  SEM.  However,  the  time-discretization  error  is  approximately 
10~n;  because  of  this,  there  is  no  gain  in  accuracy  with  further  grid  refine¬ 
ment  from  1600  to  6400  elements  unless  a  more  accurate  time  discretization 
method  is  used. 

Atmo-3D.  The  transport  of  a  passive  tracer  in  a  neutrally  stratified  atmo¬ 
sphere  in  a  large  domain  represents  an  idealized  application  to  a  seemingly 
real  atmospheric  problem.  This  final  test  is  a  proof-of-concept  to  verify  the 
behavior  of  the  methodology  over  larger  time  and  spatial  scales  that  are  of 
relevance  for  real  applications. 

The  velocity  field  is  no  longer  uniform  and  constant,  but  varies  non-linearly 
in  space  and  time  during  the  evolution  of  a  rising  thermal  perturbation  orig¬ 
inally  centered  at  the  central  lower  region  of  the  domain.  The  difficulty  of 
the  test  is  expressed  by  the  transient  character  of  the  velocity  that,  in  the 
first  instant  of  the  motion,  greatly  affects  the  stability  of  the  solution  of 
the  advection  equation.  The  problem  is  defined  as  follows  [62],  The  do¬ 
main  extends  within  Q  =  [0, 1000]  x  [0, 1000]  x  [0, 1000]  m3.  It  is  divided 
first  into  10,  and  then  into  20  spectral  elements  of  order  4  along  x  and 
z,  with  1  element  along  y.  The  simulations  final  time  is  t  =  600  s.  A 
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neutral  background  state  at  uniform  potential  temperature  8q  =  300  K  is 
perturbed  by  a  cylindrical  thermal  bubble  of  radius  rc  =  250  m,  centered  in 
(xc.  yc,  zc)  =  (500,  500,  350)  m,  and  defined  by 


O' 


A 


1  +  cos 


(34) 


where  r  =  \J(x  —  xc )2  +  (z  —  zc)2  and  A  =  0.5  K.  The  top,  bottom,  left, 
and  right  boundaries  are  modeled  as  non-penetrating  solid  walls,  while  pe¬ 
riodicity  is  imposed  on  the  front  and  back  boundaries  (y-direction) . 

The  thermal  problem  is  modeled  by  the  Euler  equations  of  inviscid  com¬ 
pressible  flows  and  solved  by  the  method  described  in  [63].  The  use  of  VMS 
to  solve  the  Euler  equations  falls  beyond  the  scope  of  this  paper,  although 
the  authors  are  currently  working  at  its  implementation  in  the  context  of 
spectral  elements.  Currently,  VMS  for  the  finite  element  solution  of  the 
compressible  Euler  equations  of  dry  nonhydrostatic  stratified  flows  can  be 
found  in  [64]. 

At  time  t  =  0  s,  the  tracer  qh  is  centered  in  the  same  position  of  6' ,  but 
is  described  by  a  cylindrical  step  function  of  maximum  intensity  Qfrma,r  =  0.5 
within  a  radius  r  <=  250  m.  The  initial  state  of  qh  and  8'  is  shown  in  Figure 
26.  After  600  s  the  rising  bubble  has  developed  into  the  structure  plotted 
in  Figure  27.  If  properly  resolved,  the  tracer  is  expected  to  have  similar 
features  given  that  the  velocity  field  derives  from  the  motion  of  the  bubble. 
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Figure  26:  Atmo-3D:  x  —  2-slice  plot  at  y  =  500m  of  the  initial  conditions  of  9'  (left),  and 
qh  (right). 
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Figure  27:  Atmo-3D:  9'  after  600  s.  Grid:  20  x  20  elements  of  order  4. 

The  velocity  field  u  is  still  until  the  warm  bubble  begins  to  move  due 
to  buoyancy.  As  soon  as  u  /  0,  the  tracer  begins  to  move  as  well.  The 
sudden  change  of  state  from  rest  to  moving  generates  oscillations  at  the 
boundaries  of  the  tracer  that  are  more  difficult  to  treat  with  respect  to  its 
analogous  steady-state  case.  Figure  28  shows  the  tracer  after  600  s  on  a 
grid  of  10  x  10  elements  of  order  4.  The  filtered  solution  and  the  solution 
obtained  with  artificial  diffusion  (Figures  28a,  28b)  have  important  under 
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and  overshoots  that  propagate  in  the  whole  domain.  The  constant  diffusivity 
coefficient  used  with  AV  is  taken  as  the  average  value  of  t  of  VMS  for  the 
same  simulation.  As  expected  from  the  previous  tests,  VMS  alone  is  not 
able  to  fully  eliminate  the  oscillations  in  the  proximity  of  the  discontinuity, 
however,  improvement  is  evident  (Figures  28c).  The  best  performance  is 
obtained  with  the  combination  VMS+DC.  Using  the  theoretical  extreme 
values  for  the  tracer  (0  <  qexact  <  0.5),  the  relative  error 


is  reported  in  Table  3.  The  same  considerations  apply  for  the  finer-grid 
solution  (20  x  20  elements  of  order  4).  The  results  are  plotted  in  Figure  29 
and  the  errors  reported  in  Table  4. 
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(b)  Artificial  diffusion  ( v  =  0.001  m2s  1) 
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Figure  28:  Atmo-3D:  Tracer  after  600  sec.  Grid:  10  x  10  elements  of  order  4. 
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Figure  29:  Atmo-3D:  Tracer  after  600  sec.  Grid:  20  x  20  elements  of  order  4. 


Table  3:  Atmo-3D:  Relative  error  e  on  the  maximum  (0.5)  and  minimum  (0.0)  theoretical 
values.  Results  for  10  x  10  elements  of  order  4. 


Method 

€min 

tmax 

AV/HV  (v  =  0.001  m2s_1) 

76.89  % 

105.37  % 

AV/HV  (y  =  0.1  m2s~1) 

61.79  % 

45.00  % 

Filter 

51.05  % 

53.53  % 

VMS 

33.67  % 

36.21  % 

VMS  +  DC 

3.53  % 

5.20  % 
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Table  4:  Atmo-3D:  As  for  Table  3,  but  for  20  x  20  elements  of  order  4. 


Method 

tmin 

£max 

AV/HV  (i/  =  0.001  m2s~1) 

65.26  % 

62.59  % 

AV/HV  (^0-lmV1) 

20.33  % 

21.06  % 

Filter 

36.38  % 

53.10  % 

VMS 

14.34  % 

11.53  % 

VMS  +  DC 

1.24  % 

4.63  % 

Remark  3:  on  the  use  of  filters  in  the  previous  results  In  the 

current  work,  filtering  was  applied  in  the  usual  way  that  has  been  used 
previously  in  SE  models  (see,  e.g.,  [44,  65,  66,  67]).  That  is,  the  filtering 
coefficients  were  defined  at  the  beginning  of  the  simulation  and  applied  af¬ 
ter  every  time-step  using  the  same  filter  matrix  for  all  elements.  It  may  be 
possible  to  obtain  better  results  with  filters  if  they  are  constructed  in  a  spe¬ 
cific  way  (e.g.,  each  element  uses  a  different  filter  matrix  that  is  constructed 
dynamically)  but  a  clear  approach  on  how  to  do  this  remains  an  open  topic 
since  this  can  be  viewed  as  a  classical  limiter  but  for  Spectral  Elements  (see, 
e.g.,  [68]). 

5.  Discussion  and  conclusions 

5.1.  Conclusions 

In  this  paper,  we  proposed  the  use  of  the  Variational  Multiscale  Sta¬ 
bilization  (VMS)  method  to  stabilize  advection-dominated  problems  solved 
with  spectral  elements.  In  the  regions  characterized  by  strong  gradients, 
we  also  combine  VMS,  a  Discontinuity  Capturing  technique  (DC),  and  the 
First-Order  Subcells  method  (FOS)  for  a  better  treatment  of  Gibbs  phe¬ 
nomena  in  the  proximity  of  boundary  and  internal  layers.  The  stabilization 
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parameter  r  that  appears  in  the  VMS  scheme  was  computed  to  include  the 
characteristics  of  high-order  spectral  elements  with  (non-equispaced)  LGL 
nodes.  Numerically,  we  demonstrated  that  this  approach  is  a  possible  al¬ 
ternative  to  the  standard  filters  used  in  the  stabilization  of  the  spectral 
element  solvers  if  the  suppression  of  unwanted  under-  and  over-shoots  is  the 
main  concern.  Stabilization  by  these  methods  is  obtained  by  introducing 
a  diffusion-like  term  that  is  controlled  and  localized  where  the  residual  is 
important  (i.e.  large  gradients).  Where  needed,  the  combined  action  of 
VMS  and  FOS  yields  encouraging  results  for  high-order  spectral  elements. 
The  algorithms  were  evaluated  on  a  set  of  standard  tests  of  increasing  dif¬ 
ficulty.  A  significant  improvement  was  observed  in  the  performance  of  the 
spectral  element  solver  as  far  as  the  control  of  extrema  is  concerned,  both 
in  the  purely  advective  and  in  the  advective-diffusive  regimes.  The  most 
important  features  of  this  new  approach  are  the  following: 

•  Unlike  hyper-viscosity,  the  subgrid-scale  diffusion  is  localized  and  con¬ 
trolled. 

•  Under-  and  over-shoots  are  greatly  suppressed  relative  to  traditional 
filters. 

•  The  VMS  method  does  not  depend  on  a  free-parameter  assigned  by  the 
user.  On  the  other  hand,  in  Algorithm  1  e  is  a  free-parameter  related 
to  the  simple  error-estimator  that  was  used.  A  more  sophisticated 
estimator  should  not  depend  on  any  user-defined  constant. 

•  Currently,  the  method  is  not  fully  mass-conservative.  This  can  be  an 
issue  for  long  term  simulations  such  as  those  for  climate  applications. 
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5.2.  Application  to  atmospheric  modeling  in  climate  and  weather  prediction 
In  [57] ,  a  Kessler  microphysics  scheme  was  implemented  within  a  spectral 
element  framework  that  requires  the  advection  of  three  moisture  variables 
(vapor,  cloud,  and  rain  mixing  ratios).  This  nricrophysics  scheme  will  be  im¬ 
plemented  in  our  Nonhydro  static  Unified  Model  for  the  Atmosphere  (NUMA) 
[63]  in  order  to  simulate  both  mesoscale  and  synoptic-scale  atmospheric  phe¬ 
nomena.  As  is  well-known,  Galerkin-based  methods  yield  1)  higher-order 
accuracy  and  2)  excellent  dispersion  properties,  which  are  both  desirable  for 
advection  schemes;  however,  the  resulting  Gibbs  oscillations  produce  strong 
gradients  that  must  be  remedied  in  some  fashion.  At  present,  a  simple- 
minded  “fixer”  is  applied  whereby  negative  values  of  the  moisture  variables 
are  set  equal  to  zero.  This  fixer  acts  as  an  effective  mass  source,  thus  violat¬ 
ing  the  conservation  properties  of  the  model.  In  addition,  this  fixer  violates 
the  function  space  that  the  spectral  element  solution  inhabits.  For  these  rea¬ 
sons,  monotonic  advection  of  tracer  variables  is  essential  for  any  atmospheric 
model.  The  proposed  VMS+DC+FOS  technique  is  a  candidate  since  it  1) 
preserves  monotonicity  better  than  the  standard  filter  approach,  2)  does 
not  significantly  increase  the  cost  of  the  spatial  discretization  scheme,  and 
3)  is  completely  local  in  nature  (i.e. ,  no  additional  communications  are  re¬ 
quired  in  a  parallel  environment),  which  is  necessary  for  scaling  on  modern 
distributed  and  hierarchical  memory  environments.  However,  the  method 
should  be  improved  for  better  mass  conservation,  especially  when  the  time 
scales  at  hand  are  large. 
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Appendix  A.  Explicit  expression  for  r  on  high-order  elements 

In  this  appendix,  we  explicitly  derive  the  expression  for  t  defined  between 
two  consecutive  LGL  points  [xigi(i),  xigi(i  +  1)].  The  bubble  obtained  from 
the  integration  of  (20)  with  boundary  conditions  b(xigi(i))  =  0  and  b(xigi(i  + 
1))  =  0  has  expression: 


x(i  +  1)  —  x(i)  ux/v  x(i)eux(l+1^u  —  x(i  +  i)eUIWA' 


x 


b(x) 


u 


The  subscript  LGL  is  omitted  to  keep  the  long  expressions  simple  to  read. 
The  evaluation  of  the  integral  (26)  yields  the  expression: 


1 


x(i)  —  x(i  +  1) 


eux(i- i-i)/i/(x(z)  —  x(i  +  l))2 


x(i  +  1)  —  x{i) 


u 


gux(i)/v  _  ^ux(i-\-\)jv 
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When  x(i) 

=  0  and  x[i  +  1)  =  h,  we  have  that 

v  h  heuh/v 

T°  ~  u2  2 u  +  euh/v  -  1  ’ 

from  which,  with  little  algebra,  expression  (23)  is  recovered: 

T  =  A  (cothfPe*)  -  pA)  • 
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